Bipolarons from long range interactions: Singlet and triplet pairs in the screened 

Hubbard-Prohlich model on the chain 
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We present details of a continuous-time quantum Monte-Carlo algorithm for the screened 
Hubbard-Frohlich bipolaron. We simulate the bipolaron in one dimension with arbitrary inter- 
action range in the presence of Coulomb repulsion, computing the effective mass, binding energy, 
total number of phonons associated with the bipolaron, mass isotope exponent and bipolaron ra- 
dius in a comprehensive survey of the parameter space. We discuss the role of the range of the 
electron-phonon interaction, demonstrating the evolution from Holstein to Frohlich bipolarons and 
we compare the properties of bipolarons with singlet and triplet pairing. Finally, we present simu- 
lations of the bipolaron dispersion. The band width of the Frohlich bipolaron is found to be broad, 
and the decrease in bandwidth as the two polarons bind into a bipolaron is found to be far less rapid 
than in the case of the Holstein interaction. The properties of bipolarons formed from long range 
electron-phonon interactions, such as light strongly bound bipolarons and intersite pairing when 
Coulomb repulsion is large, are found to be robust against screening, with qualitative differences 
between Ifolstein and screened Frohlich bipolarons found even for interactions screened within a 
single lattice site. 

PACS numbers: 71.38.Mx, 71.38.Fp, 71.10.Fd 



I. INTRODUCTION 

Electron-phonon interactions have significant effects 
in quasi one-dimensional solids and can lead to the 
formation of polarons. As shown by angle resolved 
photo-emission spectroscopy, the quasi-particles of the 
one-dimensional cuprate SrCu02^ and the conductor 
Ko.sMoOa^ are significantly modified by interactions 
with phonons. Strong electron-phonon interactions are 
also expected to be relevant in low dimensional organic 
semiconductors'^. The large phonon mediated attractions 
in low-dimensional materials may also lead to pairing of 
polarons into bipolarons. In one-dimensional (ID) sys- 
tems, bipolarons can have either singlet or triplet sym- 
metry, and the presence of triplet bipolarons has been 
observed in disordered conducting polymers using elec- 
tron spin-resonance measurements^. It is the purpose of 
this article to present a comprehensive set of simulations 
of the properties of ID bipolarons that are bound with 
long-range electron-phonon coupling. 

The polaron problem was originally formulated by 
Landau and Pekar to describe the effects of lattice po- 
larization on the properties of electrons^. As an electron 
moves through a material, it polarizes the underlying 
medium of ions and electrons leading to a dipolc field 
and thus an electron-phonon interaction. Typically the 
interaction causes a cloud of phonons to follow a single 
electron, and the combination of the particles is called a 
polaron. The mass of polarons may be significantly larger 
than that of the individual electrons^. It is now well es- 
tablished that polarons are formed in many materials, 
and are probably one of the most common excitations in 
the solid state. 

More controversial is the existence of bound pairs of 



polarons known as bipolarons. If the electron-phonon 
interaction is sufficiently attractive, then it is possi- 
ble for two polarons to pair—. It is necessary that 
phonon-mediated attraction must overcome the inherent 
Coulomb repulsion between electrons in order for the pair 
to form. Studies of the continuum limit of the bipolaron 
have shown that bipolaron may be stable, but only in 
a very small region of the parameter space of retarded 
attraction and Coulomb repulsion*. 

With the goal of understanding the properties of bipo- 
larons, we have developed a quantum Monte Carlo tech- 
nique for numerically exact simulation of two electrons 
interacting via both Coulomb repulsion and phonon me- 
diated attraction^'^. Our approach is formulated in con- 
tinuous time, and therefore has the advantage that there 
are no finite size effects from the Trotter decomposition. 
We are also able to consider general long range forms for 
the electron-phonon problem. In ID, the ground state 
properties of bipolarons forming through local and near- 
est neighbor interactions can be found using an advanced 
variational technique which predicted that long range 
electron-phonon interactions may lead to light pairs in 
j^j-jio,!! -y^g have subsequently determined that bipo- 
larons interacting via very long range attractive poten- 
tials on triangular lattices may be exceptionally light^ii^. 

Our algorithm can compute numerically exact bipo- 
laron dispersions. Previously, dispersions of the 
Hubbard-Holstein model have been computed on small 
lattices, using exact-diagonalization techniques at rela- 
tively strong coupling^^. The intermediate coupling dis- 
persions of a modified Hubbard-Holstein model have been 
computed using an approximate variational technique^"*. 
In principle, approximate dispersions could also be ex- 
tracted from the spectral functions in Ref. [isl . Bipolaron 
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FIG. 1: Example paths in ID, showing the notation and concepts of the bipolaron simulation, (a) Paths in the direct 
configuration. The ends of the paths are separated by a distance A. There are N+a kinks and N-a antikinks on the 1st path 
and equivalent numbers for the 2nd (B) path, (b) In the exchanged configuration, the paths cross. A feature of the ID problem 
is that there is always at least one shared segment between kinks in the exchanged configuration. Thus if the Coulomb repulsion 
is infinite, there can be no exchanged configuration, (c) There is a third configuration where the endpoints of the paths sit on 
the same site, which we call the ambiguous configuration. When the singlet properties are calculated, all path configurations 
are associated with a factor 1. For triplet properties, direct configurations have a sign 1, exchanged a sign -1 and ambiguous 
have sign 0. 



dispersions computed by El Shawish for the Hubbard- 
Holstein model using the same method as in Refji^ were 
also presented in Ref. [l^. 

Our aim here is to simulate the singlet and triplet pairs 
that form in systems described by a generic Hubbard- 
Frohlich Hamiltonian with local Coulomb repulsion and 
long-range electron-phonon interaction^. A crucial fea- 
ture of this model Hamiltonian is that it treats both 
strong electronic repulsion and attraction from electron- 
phonon interaction on equal footing. The Hamiltonian 
has the following form, 
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Each ion has a displacement ■ Sites labels are n or m 
for electrons and ions respectively, c operators annihi- 
late electrons. The phonons are Einstein oscillators with 
frequency to and mass M . {nn') denote pairs of nearest 
neighbor sites, and Pm are the momentum operators of 
the ions. The instantaneous interaction V{n,n') has an 
on-site component U which we always treat as repulsive 
and an off-site component V . We consider 3 forms for the 
electron ion force function, /. The Holstein interaction 
is fmifi) = KSnm, leading to a Hubbard-Holstein model. 
The near-neighbor model introduced by Bonca and Trug- 
man has the form fm{n) = K[<5„,m+i/2 + '5n,m-i/2]/2. 
The Frohlich interaction force is long-range, fmin) — 

K [(m — n)^ -I- l] '^^^ exp(— |m — n\/Rsc), where k is a 
constant and Rsc is the screening radius. It is useful to in- 
troduce the dimensionless constant A = /2Muj'^W (po- 
laron shift normalized by the half band width, W = 2t) 



which characterises the strength of the electron-phonon 
coupling. 

The Lang-Firsov transformation can be used to un- 
derstand the extreme anti-adiabatic limit of electron- 
phonon models, mapping onto a Hubbard-like system 
with instantaneous interaction between electronsi^. The 
Hubbard-Holstein model maps directly to a Hubbard 
model with a new interaction U = U2WX. Since pairs 
are bound in ID for any attractive interaction, the criti- 
cal lambda for singlet formation is Ac = U/2W. 

In the anti-adiabatic limit, the near neighbor model 
becomes a UV model with U = U — 2WX and V = 
2W^A(/)(a). The near-neighbor model may be solved 
analytically to determine the binding condition V > 
2Ut/{U + U) for singlet and F > 2t for triplet pairing^^ 
(we use the symbols U and V when discussing the UV 
model to distinguish from the U in equation [1]). Singlet 
and triplet pairs are degenerate for very large U . Since 
the wavefunction of triplet states has a node at zero sep- 
aration, triplet states do not depend on U. The high 
phonon frequency limit of the Frohlich model has a long 
range instantaneous attractive tail, which makes deter- 
mination of the binding condition quite complicated^^. 
Here, we will investigate bound states for finite phonon 
frequencies. 

A feature of the models with long-range electron- 
phonon interaction is that several different configurations 
of bipolarons can pair. Bipolarons may be of a singlet 
type, strongly bound into an onsite pair denoted SO, or 
an intersite pair separated by one lattice site denoted SI, 
which has a zero in the pair wavefunction for zero separa- 
tion. Triplets are never formed in the Hubbard-Holstein 
model, but are permitted in the models with long range 
interaction. 

This article is organized as follows: The continuous 
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time quantum Monte-Carlo algorithm for bipolarons is 
outlined in section [III In sections Hill and [TVl we simulate 
singlet and triplet bipolaron properties. Dispersions and 
densities of states are presented in section |Vl We sum- 
marize in section I VII The appendix contains complete 
details of the Monte-Carlo algorithm. 



method for simulating the polaronS. An integration over 
phonon degrees of freedom leads to an effective action. 
In the two particle case, the action is a functional of two 
polaron paths in imaginary timei^, 



II. ALGORITHM 



The continuous time quantum Monte-Carlo algorithm 
for bipolarons is an extension of a similar path-integral 



A[r{r)] = fj, (e-(^/-l-^'l) + e-(^/-l-^'l)) J] r,(r')] (1) 
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The retarded interaction mediated by the 
phonons is characterized by 'I'Aj-i^l'''), = 
J2mfm[r{T)]fm+Ar[r{T')], where the vector 
Ar = r(/3) — r(0) is an offset between the end 
configurations. Ar enables the computation of mo- 
mentum dependent properties^, and is schematically 
introduced in figure [1] The indices i and j = 1,2 
denote which fermion path is being considered. The 
dimensionless variables iD = huj/t and (3 — t/kBT. 
An instantaneous Coulomb repulsion, y[ri,r2], is also 
included. The weight of a configuration is given by 
exp(yl). 

Our electron paths are continuous in time, with hops 
between sites (kinks) introduced or removed from the 
path on each Monte-Carlo step. In contrast to the one- 
particle algorithm, kinks and antikinks must be updated 
in pairs to maintain boundary conditions in imaginary 
time. Another significant difference between one- and 
two-particle algorithms is particle exchange, leading to 
bipolarons with singlet and triplet symmetries. Measure- 
ments of the ground-state singlet bipolaron are not sub- 
ject to sign problems. It is also interesting to measure 
the properties of triplet bipolarons, which depend on the 
sign that the wavefunction picks up under exchange. Es- 
timators for the triplet bipolaron will be discussed later. 

We use a correlated weighting scheme for kink addition 
and removal, and a simple weighting scheme for path 
selection and kink type. 

Rl Choose a kink with direction (type) I from all pos- 
sible hops with equal probability Pi = 1/Nk and 
determine the anti kink as —I. Nk is the number 
of possible hops between sites (normally the num- 
ber of nearest neighbors). Pi always cancels in the 



balance equations. 

R2 Assign path A with probability 1/2 from both 
paths, and assign path B as the other path. 

R3 Choose shift type for the first kink with equal prob- 
ability Ps = 1/2. There are two shift types, the 
path above the kink insertion point can be shifted 
through Z, or the path below the insertion point can 
be shifted through —I. 

R4 When selecting the first kink to remove (which has 
type I), do so with probability 1/Nai{I)- Where / 
is the configuration before removal. That kink is at 
time T. 

R5 When selecting the second kink of type I to 
remove from path A, do so with probability 
p{T',T)/J2fJl''^^ Pin,T). Where / is the config- 
uration before removal. 

R6 Always insert the first kink at r with probability 
density p{t) = 1/(3. 

R7 Always insert the second kink at r' with probability 
p{t',t). 

R8 Choose shift type for kink 2 from kink 1, depend- 
ing if shifts are correlated (no change in interpath 
distance) or anticorrelated. 

The probability p{t' , r) can be chosen either as unity, 
recovering previous insertion and removal rules^, or 
p{t',t) can be larger for kinks that are separated by 
smaller time differences, potentially leading to improved 
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acceptance ratios at strong coupling. We use the weight- 
ing p{At) = 9{a - At) /2a + 1/2(3, where 9 is the step 
function and a is a parameter that can be made smaher 
to improve acceptance ratios (At = t' — t is mapped 
onto the interval [0,/3) and < a < /3). 

In the two particle case, it is necessary to operate on 
two kinks simultaneously to ensure that the end config- 
urations of the paths satisfy the boundary conditions. 
We have previously discussed how the algorithm works 
without exchange (applicable to ladder systems where 
particles sit on different legs). Four extra Monte-Carlo 
update rules can be used in the exchanged configura- 
tion, bringing the potential number of binary updates to 
eight. The new rules are required to ensure ergodicity in 
the exchanged configuration. Here we summarize all the 
possible updates: 

I Two kinks of the same type I are added to (or re- 
moved from) two different paths. 

II A kink-antikink pair is added to (or removed from) 
one of the two paths. 

III A kink of type I is inserted into one path, and 
another kink of the same type I is removed from 
the same path. This type of update shifts kinks in 
imaginary time. It is not essential. 

IV A kink of type I is added to one path, and an an- 
tikink —I is removed from the other path (this is 
also not essential). 

V Addition (or removal) of kink I and antikink —I on 
different paths. 

VI Kink type I is inserted on one path and kink type 
I is removed from the other. 

VII Addition of kink I and removal of antikink —I on 
the same path. 

VIII Addition (or removal) of a pair of kinks of type I 
on a single path. 

In the direct configuration, updates (I-IV) may be 
used, and the shift types of the two kinks may be identical 
(correlated) or opposite (anticorrelated) . It is convenient 
to think of the seperate correlated and anticorrelated up- 
dates as 8 different updates. In the exchanged configura- 
tion, updates (I-IV) are used unchanged, except that the 
shifts on both paths must be correlated (i.e. both top, 
or both bottom) to maintain the correct boundary con- 
ditions in time, thus I-IV only count as 4 updates. The 
remaining 4 updates (V-VIII) are used only with anticor- 
related shifts, thus in total, there are 8 possible updates 
in the exchanged configuration. Updates (V-VIII) are 
essential to access the whole configuration space as they 
change the interpath distance A of the exchanged config- 
uration. Finally in the ambiguous configuration, all the 
possible updates may be used. To simplify the scheme 
and its testing, we use a minimal update set of I, II, V 



and VIII. Details of the probabilities for binary updates 
are given in the Appendix. 

The algorithm can be made more efficient by introduc- 
ing Monte-Carlo updates to change the paths from direct 
into exchanged configurations. In the path- integral for- 
malism, an exchange of particles corresponds to swap- 
ping the ends of the paths at r = /3. An example of an 
exchanged configuration is shown in Fig. [1] To make 
the exchange in ID, the t = (3 end of path B must be 
shifted by —A lattice sites, and that of path A through 
-l-A sites. This can be achieved by inserting A — n kinks 
and removing n antikinks from path A, and inserting 
A — m antikinks and removing m kinks from path B. 
Details of how to take the continuous time limit of this 
type of update are provided in the appendix. This kind 
of exchange update is relatively slow, and we attempt 
this with a small probability of the order of 10% on any 
Monte-Carlo step. 

There is another way of carrying out exchange up- 
dates. If the paths occupy the same lattice site at the 
same point in imaginary time, then the paths have a 
segment in common. The paths can be broken at that 
segment, with the bottom of path A attached to the 
top of path B and the bottom of path B to the top of 
path A. This type of exchange can speed up the com- 
putation of triplet properties since the update is very 
fast and it can be attempted regularly without loss of 
speed. In our algorithm, we attempt a common segment 
exchange with probability pes every time an exchange is 
attempted. Once the attempt is started, we test for the 
existence of a common segment. If there are no shared 
segments, then the update is rejected and no further up- 
date is attempted on that Monte-Carlo step. Otherwise, 
the common segment exchange is accepted with proba- 
bility P{C D)= min{l,exp(A(L>) - A{C))]. 

There is another update that can be useful if the action 
has a similar value for path configurations with small and 
large separation (as is typically the case when the bipo- 
laron is only just bound). Then it can be appropriate to 
attempt to shift one of the paths through a longer dis- 
tance when the paths are in a direct configuration. This 
update enables more rapid sampling of the path configu- 
rations and is especially useful for making accurate mea- 
surements of the inter-particle separation. 

After a warmup period, we make a series of measure- 
ments. The measurements are made every few updates 
to avoid correlation between Monte-Carlo steps. We also 
use a blocking procedure in conjunction with bootstrap 
resampling for error estimation. In situations where the 
measurements take a long time (comparable with the 
computation of the action) - such as the estimators for 
the total energy and isotope exponent, this sparse sam- 
pling acts to speed up the algorithm, since correlated 
measurements do not improve estimation of the averages. 
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Estimators relevant to the current manuscript are the 
total singlet energy, 



E 



lim 




(2) 



where Ni is the number of kinks of type i, and angu- 
lar brackets denote ensemble averaging. The number of 
phonons associated with the bipolaron, 



^ph = - lim 



1 



dA 



I3{s) \ duj 



(3) 



where the derivative is taken keeping \u) constant, s ~ 
1 when considering singlet states, and s = 1,0,-1 for 
measurements of triplet states, as described in Fig. [T] 
The bipolaron band energy spectrum can be computed 
from 



e(fc) - e(0) = - lim -An 



[s cos(fc • Ar)) 



(4) 



where k is the quasi momentum. By expanding this ex- 
pression in small k, the i-th component of the inverse 
effective mass is found to be 



1 



— lim 

TO* /3^oo 
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(5) 



We compute the bipolaron radius as the root mean square 
distance between paths. 



R. 
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(r(T)i-r(T)2)2dr) (6) 



The mass isotope coefficient, a„i» — dhim* /dlnM , is 
calculated as follows 



lim — - 
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(8) 



The singlet-triplet splitting can also be determined, 

A.* = -^ln((.)) 

Since (s) is always less than 1, the singlet-triplet split- 
ting is always positive and the energy of the triplet state 
is higher in energy, as expected since the ground state 
wavefunction of two particles should have no nodes. 

An alternative way of computing spectra and density 
of states is also available. We take the measurement 
{dr,Ar) for values of r consistent with a few lattice spac- 
ings. The magnitude of the measurement drops off very 
rapidly with r. Then the whole spectrum can be com- 
puted for any k point, without deciding which k to in- 
vestigate before the algorithm starts. The expression for 
the dispersion calculated in this way is. 



e/e - So = - 111 



' cos(fc.r) 



\SSr,Ar) 

is) 
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This form is particularly useful for determining the den- 
sity of bipolaron states. 

In the following, we have restricted the electron paths 
to remain withing 200 lattice sites of each other by en- 
forcing a pairwise infinite instantaneous potential well 
dependent on the inter-particle separation. This stops 
paths from moving apart indefinitely. The 200 site well is 
sufficiently large to ensure that finite size effects are neg- 
ligible compared with other considerations such as the 
statistical error in the Monte Carlo averaging. 



III. SINGLET BIPOLARON 

In this section, we examine and compare the properties 
of the singlet bipolarons formed in the Hubbard-Holstein 
model, screened and unscreened Hubbard-Frohlich mod- 
els and in the nearest neighbor model. We compute the 
total energy, inverse mass, inverse radius, mass isotope 
exponent and number of phonons associated with the 
bipolaron cloud. 

We begin by examining the total energy associated 
with the singlet bipolaron, which can be seen in Fig. [21 
The energy is plotted for a range of U/t and A with fixed 
Lj = 1. In the Hubbard-Holstein model, there is a rapid 
change in the gradient of the plot around the critical cou- 
pling for binding (for large U, the graph is flat because 
the two polarons are not bound, and so do not interact 
through the local U). The long range tails in the screened 
Hubbard-Frohlich model with Rgc = 3 lead to energies 
that are similar to the near-neighbor model and the dif- 
ference in energies is barely detectable by eye. This is 
a strong indicator that the nearest neighbor interaction 
strength ^(a) is the most important factor in determin- 
ing the propoerties of the bipolaron. ^(a) = 0.5 in the 
near neighbor model and = 0.4688 to 4 significant 

figures in the Rsc = 3 screened Frohlich model (a is the 
lattice vector) . The unscreened Frohlich interaction leads 
to bipolarons that are more strongly bound. For exam- 
ple, a much larger U/t is needed to unbind the bipolaron 
when the Frohlich interaction is simulated at A = 0.3. In 
the limit of low screening radius, the energy of the bipo- 
laron approaches the Holstein limit very slowly, with the 
very rapid binding close to a critical U associated with 
the Holstein model not seen in models with long range 
interaction. This is a qualitative difference between Hol- 
stein models and those with a small amount of intersite 
attraction. 

In the very strongly bound limit, the effects of electron 
hopping are small. An effective atomic Hamiltonian can 
be found by applying the Lang-Firsov transformation to 
the phonon part of the Hamiltonian, 
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Reintroducing the effects of Coulomb repulsion, the total 
energy of the strongly bound onsite bipolaron (formed 
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FIG. 2: Total energy associated with the singlet bipolaron for the 3 different models. The limits corresponding to strongly 
bound SO bipolarons are shown as dotted lines, and arrows on the right axis correspond to strongly bound SI bipolarons. The 
surprising result here is that the long range tails in the Hubbard-FrohUch model with Rsc = 3 lead to similar results to the 
near-neighbor model proposed by Bonca and Trugman. In the results for the Hubbard-Holstein model, there is a sudden change 
in the gradient of the plot around the critical coupling. This is not apparent for either of the models with long range coupling. 
/3 = 14 in all figures in this section. 



when 8tX > U) can be computed to be, 

E = U-8tX (11) 
The energy of the strongly bound SI bipolaron at large 



A and U > 2tX is. 



£;= -4tA[l + $(a)/$o]. (12) 



These limits are also shown in Fig. [5] (the strong bind- 
ing SI energies are only shown for A = 0.9, 1.2 and 1.5). 
There is a clear convergence between the numerics and 
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FIG. 3: Total number of phonons associated with the singlet bipolaron. The limits corresponding to strongly bound SO 
bipolarons are shown as arrows on the left axis, and arrows on the right axis correspond to strongly bound SI bipolarons. The 
similarity between the two models with long range interaction and the distinction with the Holstein bipolaron can immediately 
be seen. The nearest neighbor and lattice-Frohlich models show only quantitative differences. In particular, the crossover from 
the weakly to strongly bound states is gentle in comparison with the Holstein interaction, where the bipolaron is suddenly 
bound at a critical coupling, causing a sudden increase in the number of phonons associated with the bipolaron. 



the strongly bound SO state at low U. Bipolarons with 
long-range interactions converge on the limit more slowly 
than Holstein bipolarons. The SI energies at large U be- 
gin to agree at large A. This is expected, since the bipo- 
larons formed in models with long-range interactions are 
quite mobile, so the efllects of hopping should be con- 



sidered, even when A — 1.5. Note that there is no SI 
bipolaron in the very large U ^ 2tX limit of the Holstein 
model. 

The total number of phonons associated with the sin- 
glet bipolaron is shown in Fig. [31 The similarity between 
the pairs in the models with long range interaction and 
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FIG. 4: Inverse mass of the singlet bipolaron. It is clear that both the more complicated lattice-Frohlich bipolaron and 
the simplified model with nearest neighbor interaction generate bipolarons with similar light mass. The inverse mass of the 
nearest-neighbor model has a small hump around U — 3. 



the distinction with the Holstein bipolaron is clearly visi- 
ble. The near-neighbor and lattice-Frohlich models show 
only quantitative differences. In particular, the crossover 
between weakly and strongly bound cases in the models 
with long-range interaction is gentle in comparison with 
that seen in the Hubbard-Holstein model; in the case of 
the Holstein interaction the bipolaron is rapidly bound 
on decreasing U, causing an abrupt increase in the num- 



ber of phonons associated with the bipolaron. 

The total number of phonons associated with the 
strongly bound onsite bipolaron (following the argument 
in Ref. |^ is, 

8A , , 

A^ph = — . (13) 

id 

Note that the number of phonons in the strongly bound 
onsite bipolaron does not depend on U, so the number of 
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FIG. 5: Inverse pair size associated with the singlet bipolaron. The Holstein bipolaron always breaks up at large U, indicated 
by the vanishing inverse radius. In contrast, all the models with long-range interaction have bound SI states at large U when 
there is sufficient electron-phonon coupling. 



phonons reaches a limiting value on decreasing U. This 
limit can be seen in Fig. [3] as the arrows on the left- 
hand y-axis. The Holstein bipolaron rapidly approaches 
the strongly bound SO limit on decreasing U, with the 
SO limit approached less rapidly as Rsc increases. The 
reason for this is that the longer range tails lead to a 
shallower effective potential, so SO bipolarons are harder 
to bind into purely on-site pairs. The number of phonons 



associated with the SI bipolaron can also be found, 

A^ph = [l + $(a)/$o]. (14) 

UJ 

This is represented as arrows on the right hand side of the 
graphs. Again it is clear that hopping effects should be 
included at intermediate A, with the agreement becoming 
better as A is increased and the bipolaron becomes more 
strongly bound. 
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FIG. 6: Mass isotope exponent of the singlet bipolaron. The dip in the exponent seen in the models with long range interaction 
corresponds to the crossover between SO to SI bipolarons. In general, the mass isotope exponent is much larger for SO pairs 
(low U) than SI pairs (high U, A and long range interactions). 



We show the inverse mass of the singlet bipolaron in 
figure m Both the more compUcated Hubbard- Frohlich 
bipolaron and the simplified model with only nearest 
neighbor interaction generate bipolarons with similar 
light mass. At large U, the mass of the Holstein bipolaron 
does not change, because the bipolaron is not bound. 
This is also true of the long range models at small enough 
A. There is a notable feature in the mass at intermedi- 



ate [/ 3 and A = 1.5. The mass associated with the 
nearest-neighbor bipolaron has a hump, with the bipo- 
laron becoming lighter on increasing U, before becoming 
slightly heavier again. 

We have not been able to compute analytical expres- 
sions for the mass in the strong coupling limit of the 
Hubbard-Frohlich model since each electron hop excites 
phonons at an infinite number of sites, complicating the 
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second order perturbation theory. For the near-neighbor 
and Holstein models, the strong binding hmits have been 
discussed by Bonca and Trugma n^"'^^ . As noted by 
Bonca and Trugman, the energies of SO and SI config- 
urations in the near-neighbor model with high phonon 
frequency are degenerate when U — W\ at large uj^. 
This also occurs in the Hubbard-Frohlich model when 
U - 2WX = 2VFA0(a). At this point, the motion of the 
bipolarons is 1st order in t since the electron does not 
need to tunnel through a high energy barrier to move, 
leading to the mass decrease that can be seen in Fig. 
2) In this sense, the mass decrease is similar in origin 
to the superlight small bipolarons seen on lattices with a 
triangular component^. The decrease in mass at interme- 
diate U is less pronounced in the case of the unscreened 
Frohlich interaction, and not visible to the eye on the 
plot, but is present, with a peak at around U/t — 4 for 
the A = 1.5 plot. 

The pair size associated with the singlet bipolaron 
(shown in Fig. [5]) demonstrates the crossover between 
SO and SI bipolarons. There is a qualitative difference 
between models with short and long range interaction. 
For larger electron-phonon couplings and long range in- 
teractions, the bipolaron size changes from small (onsite) 
Rbp < a to intersite Rbp ~ a (the lattice constant) on 
increasing U. The change between these two types of 
bipolaron occurs at [/ ~ 3 for A = 1.5 in the nearest- 
neighbor model. The pair size associated with the un- 
screened Frohlich interaction is not as small as that for 
the screened Frohlich interaction, presumably because 
the potential well does not change so rapidly with lattice 
index. It can be seen that at large U, the SI bipolaron be- 
comes more strongly bound on increasing A, with the in- 
verse radius tending towards unity a.tU— 14. Again this 
demonstrates that the SI bipolaron can only be consid- 
ered to be strongly bound at large A and U ^ X. In con- 
trast, the Holstein bipolaron is always unbound at large 
[/, as the bipolaron radius becomes infinite {R^p 0). 
The modification of the binding diagram on changing Rsc 
can be seen in these figures. For Rsc = 1, A = 1.5 is 
needed to bind the bipolaron at strong U, whereas the 
Frohlich bipolaron binds at A = 0.9. 

We finish this section by computing the mass isotope 
exponent, a, which is plotted in Fig. [5] for a variety of A 
and U. A dip in the isotope exponent seen in the nearest 
neighbor and screened Frohlich interactions corresponds 
to the hump seen in the inverse mass. Again, this dip 
coincides with the degeneracy in the SO and SI energies. 

For completeness, we note that it is also possible to 
compute properties of the bipolaron at weak coupling 
and U > where neither singlet nor triplet bipolarons 
are bound, i.e. the measured properties correspond to 
those of two polarons. The weak coupling properties of 
the polaron can be determined from second order pertur- 
bation theory and are found to be E/t — —z — XrE{(^), 
mo/™* = 1 — Arm*(w), = AFfj^^. ((!>) and Np^^ — 
^^Npt,i<^)- For (D = 1, the F coefficients can be found 
in Ref. [2l|, with the exception of the coefficients for 
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FIG. 7: Top: Total triplet and singlet energies of Hubbard- 
Frohlich bipolarons. U = 14t, /3 = 7, w = 2 and U/t = 14, 
uj/t = 4, (3 — 7. Bottom; Singlet-triplet splitting when (3 — 7, 
ijj = 2 and various U/t and A. 



the nearest neighbor model which are F^; = 2.391(5), 
r,„. = 0.607(7), r^^, = 0.294(1) and F^Vp,, = 1.999(1) 
and for the Hubbard-Frohlich model with Rgc = 2, which 
are F^ = 1.423(7), F„. = 0.392(0), F^^. = 0.1540(0) 
and Fat^j, = 1.116(2). 



IV. TRIPLET BIPOLARONS IN THE 
HUBBARD-FROHLICH MODEL 

In this section, we compute the properties of triplet 
bipolarons and compare with those of singlet bipolarons. 
Our algorithm can compute triplet and singlet properties 
simulaneously since they are related through the mea- 
surement of the exchange sign. Strictly, it is not neces- 
sary to directly calculate the properties of triplet states 
in ID in this way (singlet and triplet bipolarons are de- 
generate at large U) however this contributes part of the 
development of our algorithm, since there is no alterna- 
tive approach for 2D and above. The Hubbard-Holstein 
model in ID has no bound triplet states, so the proper- 
ties associated with triplet symmetries are not examined 
for local interactions. 

A quantity that is straightforward to measure using our 
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FIG. 8; Top: Inverse mass of singlet and triplet bipolarons 
when /? = 7, il) = 2, [/ = 14t. The Hubbard-Frohlich model is 
simulated. N.B. Triplets are heavier than singlets at strong 
coupling, as in the solution of the UV model. Bottom: Vari- 
ation of inverse mass with U . The triplet mass is constant. 




FIG. 9: Number of phonons associated with singlet and triplet 
bipolarons in the Hubbard-Frohlich model. Inset: the differ- 
ence between the number of phonons associated with singlet 
and triplet bipolarons. These data were computed for the 
paramters /? = 7, a; = 2 and V = 14t. 
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FIG. 10: Inverse size of singlet and triplet bipolarons in the 
Hubbard-Frohlich model when j3 ^ 7, uj = 2, U = lU 



Monte-Carlo algorithm is the energy splitting between 
singlet and triplet bipolaron states. From this and the 
singlet energy, the total energy of triplet bipolarons can 
be determined. The singlet-triplet splitting and singlet 
and triplet energies are shown in figure [T] The top panel 
shows the total singlet and triplet energies for various 
A and Cu — 2 and U — 14. Below a certain electron- 
phonon coupling, the singlets and triplets are degener- 
ate, as Coulomb repulsion inhibits pairing in the singlet, 
and the triplet state (which must be higher in energy 
than the singlet) is unbound. This can also be seen for 
uj/t — 4, and U/t ~ 14. The singlet-triplet splitting 
is relatively small compared to the total energy of the 
bipolaron and the total energy only changes slightly on 
increase of phonon frequency. 

The bottom panel of Fig. [7]shows the splitting between 
singlet and triplet states for a variety oiU/t and A. On 
decreasing U , the singlet-triplet splitting increases. This 
is expected as the singlet state becomes increasingly well 
bound on decreasing U , whereas the binding energy of the 
triplet should remain constant. Efficient computation of 
the splitting is limited by the temperature, as the higher 
energy triplet configurations need to be visited enough to 
obtain accurate statistics. For the (3 = 7 shown here, this 
effectively limits us to maximum splittings of the order of 
t/2 (otherwise computation time is excessive). Increase 
of A also increases the splitting, which can also be seen 
in the top panel. 

We also examine the relative masses of the triplet and 
singlet bipolarons, which can be seen in Fig. [51 We be- 
gin by computing the inverse effective mass for w = 2, 
U/t = 14 and (3 = 7 which can be seen in the top panel 
of the figure. For large A, it is clear that the triplet bipo- 
laron is heavier than the singlet one. However, for small 
A, the situation is reversed, with the triplet state slightly 
lighter than the singlet one. The slightly lighter mass 
is persistent to A ~ 1.125, where the masses cross, for 
large A, triplets are observed to be heavier than singlets. 
We also make computations for U/t = 14, and (D = 4. 
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Triplets are heavier than singlets over a wider range of 
the parameter space for larger lo. This is consistent with 
the large lo limit, where the models with long range in- 
teraction have similar properties to the UV model (the 
UV model always has heavier triplets). 

The bottom panel of Fig. [5] shows how the masses of 
singlets and triplets change with U . On decreasing from 
large [/, the singlet mass decreases to become lighter than 
the triplet as the SO and SI states become degenerate. 
However at small U this is reversed, since the singlet 
becomes tightly bound in an SO state (see Fig. 2). This 
binding into and SO state, and the associated increase in 
mass can also be seen in the A = 0.5 curves at low U . 

The number of phonons associated with singlet and 
triplet bipolarons is shown in figurelHl with the difference 
between them shown in the inset. Simulations were run 
ioi P — 1 , u} — 2, U = 14i and various A. The number of 
phonons associated with singlet and triplet bipolarons is 
very similar. For the parameters shown, there are slightly 
more phonons associated with the singlet than the triplet, 
but the number of phonons is degenerate when A is very 
small (as the two particles are not bound). The number 
of phonons associated with the cloud is also similar for 
singlet and triplet bipolarons at very large A. For such 
a large U and A = 2, the singlet bipolaron has some SI 
characteristics which explains the similarity between sin- 
glet and triplet. For much larger X:^ U/t, the electron- 
phonon coupling is expected to overcome the Coulomb 
repulsion to bind the singlet bipolaron into an SO con- 
figuration, and the properties of singlet and triplet will 
differ. 

We complete this section by considering the inverse size 
of singlet and triplet bipolarons which can be seen in Fig. 
[TUl As before, /3 = 7, <Z> = 2, and C/ = Ut. The triplet 
is larger for all A. This makes the heavier triplet states 
slightly counterintuitive, because in the case of singlets, 
larger wavefunctions tend to be associated with lighter 
bipolarons. 



V. BIPOLARON DISPERSION 

One of the powers of our Monte-Carlo algorithm is 
the possibility to compute the dispersions of the bipo- 
laron efficiently. In this section, we discuss spectra com- 
puted with our algorithm. The estimator for our spectra 
is given in equation [4] We compute triplet and singlet 
spectra at the same time. 



A. Hubbard— Holstein model 

In the Hubbard-Holstein model, triplets are never 
bound, so we do not compute the dispersion associated 
with triplet symmetry, concentrating instead on singlets. 
We show the effects of varying U/t on the dispersion of 
the Hubbard-Holstein bipolaron for A = 1, (D = 1 in Fig. 
1111 On decreasing U, the electron-phonon interaction 
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FIG. 11: Top: Dispersion of the Hubbard-Holstein for A = 1, 
Lj = 1, = 14 and various U. Bottom: Comparison between 
the shape of the dispersion, e(fe)/e(7r) and the bare cosine 
band. 
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FIG. 12: Dispersion of the Hubbard-Holstein model for A = 2, 
U/t = 5, /3 = 7 and various ui. Errorbars are too small to be 
visible. 
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FIG. 13: Dispersion of the Hubbard- Holstein model for cu/t ■ 
4, U/t = 5, P/t = 3.5 and various A. 



begins to dominate, resulting in strongly bound on-site 
pairs. This is characterised by a significant decrease in 
the band-width, which is consistent with a large increase 
in the effective mass. For U/t }^ 5, the dispersion does 
not change significantly on increasing U (the bipolaron is 
not bound at larger U). This is consistent with the esti- 
mate from the large frequency limit of the model, where 
bipolarons unbind if U > AtX. 

The bottom panel of Fig. [TT] shows the shape of the 
dispersion, e(fc)/e(7r) as compared to the cosine disper- 
sion associated with the non-interacting limit. The dis- 
persion is skewed away from the cosine and is flatter at 
larger momenta. The dispersion also becomes shallower 
at larger U where the bipolaron unbinds into two po- 
larons, demonstrating that the band flattening in bipo- 
larons is significantly larger than in polarons. 

The effects of varying uj are shown in Fig. [T^l As the 
phonon frequency is increased, the bipolaron band gets 
wider, which is consistent with the decrease in effective 
mass. The change in band width is quite significant as 
w/t changes from 2 to 4, and is around a factor of 10. 
We also compute the dispersion of the Hubbard-Holstein 
model for cj/t = 4, U/t ^ 20 and various A, which can 
be seen in Fig. [THl On increased A, the band width 
decreases, consistent with the significant increase in ef- 
fective mass. The bipolaron binds into an SO state on 
increasing A, with the bandwidth changing by more than 
three orders of magnitude from A = 1 (not shown in fig- 
ure) to A = 4. We note that similar calculations were 
presented in Ref. [1^. 



B. Hubbard— Frohlich model 

The form of bipolaron dispersions in the Hubbard- 
Frohlich model has not previously been computed. At 
large coupling and large phonon frequency, the Hubbard- 
Frohlich model maps onto a UV-model, which is a conve- 
nient starting point for analytics. The singlet dispersion 



of the UV model can be found from the equationi^, 

(Ves(fc)2 - 16^2 cos2(fc/2) + L/) X (15) 
(Ve3(fc)2 - 16*2 cos2(fc/2) - es(fc)) - 2V {U - e,(fc)) = 



the triplet dispersion is given by 



4*2 



et{k) ^ -V-^ cos" {k/2). 



(16) 



here, we take U = U — 2W\ and V = 2WX(j){a), which is 
an acceptable approximation for a very tightly bound 
state of the Hubbard-Frohlich model at very large u; 
(since at very large A the magnitude of the pair wave- 
function tends to zero too rapidly to sample next-nearest 
neighbor interactions). 

We plot dispersions for various oj/t with U/t — 20, 
A — 8, (3/t — 3.5 and a Frohlich interaction in figure 
[131 All dispersions are plotted relative to the zone center 
singlet energy. We also plot the UV model dispersions 
on the panel showing the uj/t = 40 curves. We have al- 
ready determined that triplet bipolarons are significantly 
heavier than singlets in the large U large A regime, and 
this is consistent with the much flatter triplet disper- 
sion. In analytics of the UV modeU^, singlet and triplet 
dispersions are found to be degenerate a,t k ~ ir/a. In 
agreement with results from the UV-mode\, the spectra 
are nearly degenerate at the zone edge for large phonon 
frequencies, but become increasingly separated in energy 
for lower phonon frequencies. This is clear evidence that 
the degeneracy at the tt point in the UV-mode\ is lifted 
by dynamical effects. 

We also investigate the effect of varying A on the sin- 
glet and triplet dispersions in Fig. [121 Data are plot- 
ted for the parameters cu/t = 4, P/t = 3.5, U/t — 20 
and for various A. Where no errorbars are visible, they 
are smaller than the points. The dispersions are plot- 
ted relative to the zone center singlet energy. The triplet 
spectrum is much flatter than the singlet one, which is 
consistent with the larger triplet masses at higher uj and 
A. In particular, the triplet dispersion is almost flat for 
A = 12 in comparison with the singlet dispersion. For 
A = 1, where singlets and triplets are only just bound, 
the dispersions are near degnenerate. In contrast to the 
Hubbard-Holstein bipolaron, the bandwidth decreases by 
only one order of magnitude on increasing A from 1 to 
4.5, which is a feature of the long-range interaction. 

We are also interested in how the dispersion varies as 
the Coulomb repulsion is increased. Variation of the sin- 
glet spectrum with varying U is shown in Fig. 1161 Data 
are plotted for parameters oj/t = 4, A = 2.5, (3 = 3.5 and 
a Frohlich interaction. The triplet dispersion does not 
vary with U , and as U is increased to very large values 
singlet and triplet dispersions converge. 

In figure [T7l we compute the zone center and zone edge 
singlet-triplet splitting as uj/t is varied. X = 8, U/t = 20, 
/? = 3.5. As phonon frequency increases, the zone edge 
singlet and triplet states become closer in energy with 
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Jt 



k [a-^] k [a-^] 

FIG. 14: Singlet and triplet dispersions. Data are shown for various w/i, U /t = 20, A = 8, P/t — 3.5 and a Frohlich interaction. 
Error bars are smaller than the points. The dispersions axe plotted relative to the zone center singlet energy. Note that 
the triplet spectrum is much flatter. The spectra are near degenerate at the zone edge for large phonon frequencies, but axe 
separated in energy for lower phonon frequencies indicating that the separation is a dynamical effect. The dispersions of the 
correponding UV model are also shown on the w/t = 40 plot. 

singlet and triplet expected to become degenerate at very bipolaron. 
large u as is found in the UV model. Contrary to this, 
the zone center singlet-triplet splitting increases with in- 
creased w. This is mainly due to decrease in singlet en- 
ergy. The triplet also becomes lower in energy on increas- 
ing u), but the decrease is smaller than for the singlet 
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FIG. 15; Singlet and triplet dispersions. Data arc plotted for the parameters, cu/t = 4, U/t = 20, various A, j3/t = 3.5 and a 
Frohlich interaction. Where no orrorbars arc visible, they are smaller than the points. The dispersions are plotted relative to 
the zone center singlet energy. Note that the triplet dispersion is much flatter than the singlet one. 



VI. SUMMARY 

In this article, we presented the results of simulations 
of bipolaron properties for a scries of clcctron-phonon 
models in ID. The models vary in terms of the inter- 
action range, from the simple on-site Holstein interac- 
tion, through the model suggested by Bonca and Trug- 
man with interaction only between particles on nearest 



neighbor sites to the full lattice Hubbard-Frohlich model. 
The continuous time quantum Monte Carlo algorithm 
was used to simulate the models, and full details of our al- 
gorithm have been presented. A wide range of bipolaron 
properties have been calculated, with a detailed overview 
of the parameter space and comparison made between 
models with differing ranges of electron-phonon interac- 
tion. In general, only a small long range interaction is 
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FIG. 16: Variation of singlet spectrum with varying U. Data 
are plotted for parameters uj/t = 4, A = 2.5, /3 = 3.5 and a 
EYohlich interaction. Errorbars are smaller than the points. 
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necessary to significantly modify model properties from 
the strictly local Holstein case, with light masses and 
strong pairing present for even relatively strong screen- 
ing. 

The strength of the nearest neighbor effective interac- 
tion is found to be the most important factor that influ- 
ences bipolaron properties. This increases the possibility 
of engineering a bipolaron condensate in higher dimen- 
sions, since with shorter tails, clustering corresponding to 
phase separation is less likely. Given that the simplified 
near neighbor model is much better suited to analytics 
than the lattice Frohlich approach, it is suggested that it 
should be adopted as a more general model to investigate 
electron-phonon problems in higher dimensions. Mod- 
els involving ions oscillating above the plane or chain of 
electrons have been suggested to derive this and similar 
forms. This type of configuration is not possible in 3D. 
However, in D > 2 a long range component to electron- 
phonon models is expected, since a perfectly momentum 
independent interaction is unlikely. Simulation of such 
models will form part of future publications. 
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APPENDIX 

The formal derivation of the continuous time Hmit for 
updates (V) to (VIII) by considering detailed balance is 
very similar to that for updates (I) to (IV), which has 
been described in detail elsewhere^^ (although here we 
make a slight modification for weighted insertions, where 
kinks are more likely to have similar values of r accord- 
ing to the weighting p(r, r')). Thus, we simply summa- 
rize the rules needed for implementation of the algorithm 
here. We use only updates I, II, V and VIII as a mini- 
mal update set that allows the algorithm to efficiently 
explore all path configurations. A detailed derivation 
is then given for the multiple kink update relating ex- 
changed and direct configurations in subsection [51 since 
there are some subtle considerations. 



1. Update (I): Kink addition on both paths 

Consider two path configurations C and D, where D 
has two additional kinks, I on path A and I on path B. 
The update proceeds by selecting kink shifts, paths and 
kinks to add and remove according to the scheme, 

(i) Choose kink types shifts and paths according to the 
probability rules R1-R3. 

(ii) Propose pair removal with probability 1/2. Other- 
wise attempt pair addition. 

(iii) If Nai = or Nbi = in the initial configuration 
and pair removal was proposed, then the update is 
rejected. 

(iv) For pair addition, select times with equal probabil- 
ity density according to rules RS] and FJS] 

(v) For pair removal, select an I kink from path A and 
a I kink from path B according to rules Fi[6]and Fi[7l 



(i) Choose kink types shifts and the path according to 
the probability rules R1-R3 

(ii) Propose pair removal with probability 1/2. Other- 
wise attempt pair addition. 

(iii) If Nai = or Na-i = in the initial configuration 
and removal was proposed, then reject the update. 

(iv) For pair addition, select times with probability den- 
sity according to rules Fill] and Fil5l 

(v) For pair removal, select an I kink from path A and 
a —I kink from path A according to rules ElHl and 
Id 



P(addition) = min < 1, 



^2^gA(D)-A(C) 



NiAiD)j: 



(A.2) 



3. Update (V): Addition of kink and antikink on 
different paths 

Consider two path configurations C and D that differ 
only in that D has an additional hop I on path A and —I 
on path B. 

Using the following rules: 

(i) Choose kink types shifts and paths according to the 
probability rules R1-R3 

(ii) Propose pair removal with probability 1/2. Other- 
wise attempt pair addition. 

(iii) If removal is chosen and Nai = or Nb-i = in 
the initial configuration then removal is aborted as 
a rejection. 

(iv) For pair addition, select times with equal probabil- 
ity density according to rules RSI and El5l 

(v) For pair removal, select an I kink from path A and 
a —I kink from path B according to rule FjEl 



P(addition) — min < 1, 



^2p^A{D)^A{C) 



Nia{D)Y: 



Nib{D) 



(A.l) 



P{C D)= min <^ 1, 



^2f^^A{D)-A(C) 



Nia{D)Y. 



'AD) 



(A.3) 



2. Update (II): Addition of kink and antikink pair 
on one path 

Consider two path configurations C and D that differ 
only in that D has additional kinks I on path A and —I 
on path A. 

Using the following rules: 



4. Update (VIII): Addition/removal of kink pair to 
a single path 

Consider two configurations C and D, where configura- 
tion D has two more kinks on path A than configuration 
C and the same number of kinks on path B. There is a 
slight complication here because there are two ways of 
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inserting kinks to reach the end configuration, so we al- 
ways use equal time probabilities on this update (i.e. no 
weighted kink insertions). 

(i) Choose kink types shifts and paths according to the 
probability rules Fl[T]-r@ 

(ii) Propose pair removal with probability 1/2. Other- 
wise attempt pair addition. 

(iii) If Nai < 2 and removal has been chosen, the update 
is aborted. 

(iv) For pair addition, select times with equal probabil- 
ity density according to rule Fi|4l We do not corre- 
late kink times on this update due to ambiguities 
with kink insertion ordering. 

(v) For pair removal, select two different I kinks from 
path A. 



If the Metropolis scheme is used, 

(A.4) 



5. Exchange by kink insertion and removal 

In the path- integral formalism, an exchange of particles 
corresponds to swapping the ends of the paths at t — [3. 
An example of an exchanged configuration is shown in 
Fig. [1] To make the exchange, the t = (3 end of path B 
is shifted by —A lattice sites, and that of path A through 
+A sites. From here on, we refer to the t = f3 end of a 
path as the 'top' of the path. This section explains how to 
take the continuous time limit of this more complicated 
class of update. 



In the conventional way, update probabilities are determined from the detailed balance equation, 

W{C) ■ g(forward) • P{C ^ D) = W{D) ■ g(inverse) • P{D C), (A.5) 

where W{I) represents the statistical weight of configuration / and Q is the probability of selecting a specific update 
from all the available updates (for example a specific hop direction is chosen and a value for the number of antikinks 
to remove). Configuration D is the exchanged configuration, and C represents the non-exchanged configuration. 

The total difference in hopping quanta between initial and final configurations is A — 2m since n — A — m kinks 
are added and m antikinks are removed. Therefore, the ratio of statistical weights in the partition function is, 

W{D) _ ,^^,A~2m A a{D)~Aa{C) 



W{C) 



^^^^■^A^2m^AMU)^AMV) (A.6) 



where the subscript A indicates that only path A is considered for now. 

One way of making the -l-A shift at the top of path A is to insert A — m kinks and to remove m antikinks. The 
maximum number of antikinks that can be removed is m = Na-i if A > Na-i (since it is not possible to remove 
more antikinks than exist) or m = A otherwise. The minimum number of antikinks that can be removed is m = 0. 
Therefore, there are Np possible choices of m, where Np = A -I- 1 if A < Na-i, otherwise Np = Na-i + 1- It is 
convenient to write the number of possible updates as the single expression, 

7Vp = min(7VA-i,A) + l. (A.7) 

The value of m is chosen with equal weighting 1/Np. 

Now that a value of m has been chosen, it is necessary to determine the probability that specific kinks are inserted 
into or removed from a specific configuration. The probability of removing m of the available Na-i antikinks is simply 
I/at^ jCto. The probability of choosing to introduce n = A — m kinks into configuration C with equal weighting 
1//3 is (A — m)!(AT//3)'^^™ (the probability is increased because there are (A — to)! possible ways to order the kink 
insertion that are not distinguishable in the final configuration). 

All of these considerations contribute to the probability of updating from configurations C to D via the product, 

(A-m)! /Ar\^-™ 1 



Q(forward) = — — (A.8) 

The inverse process changes the path from configuration D to C . Similar to the forward process, up to A kinks 
can be removed if they exist, otherwise, only Nai -I- A — to = Nai + n kinks may be removed. Therefore there are 
Np = m\Ti[NAi -l- A) -f 1 possible ways to choose how many kinks to remove. In the reciprocal process, to antikinks 
are inserted, and A — m kinks are removed from Nai + A — m kinks, so the total probability for choosing a specific 
inverse process is, 

/At\" 1 

Q(mverse) = — ————— — (A.9) 

m\n[N Al + n, IX) + I \ P J NAi+A-mCA-m 



20 



Substituting these probabilities into the detailed balance equation, it is possible to write the ratio update probabilities 
for a single path as, 

PAiC^D) _ m! min(jV^_^A) + l ( ^A'""'"^ .A-2m n.-.C^ . a in^ 

Pa{D -> C) (A - m)! min(iV^j +n,A) + l\ p J ^ ' N^,+A-mCA-m ^ ' ' 



There is a cancellation of the discrete time interval At, 
so the continuous time limit can be taken immediately, 

Pa{C ^ D) ^ (.g.2n-A ram{NA-i, A) + I n^.^Pa-u 
Pa{D^C) min(7V^j + n. A) + 1 Ar^,+„P„ 

(A.ll) 

A similar result is obtained for path B, with kinks re- 
placed by antikinks and vice-versa. The probabilities for 
the two paths are independent, and so the total proba- 
bility is simply a product of the two. Once the ratio of 
update probabilities has been calculated, the Metropo- 
lis update scheme can be applied to obtain a suitable 
Monte-Carlo procedure. We note that the factorials in 
this expression are a potential cause of overflow, so we 
work with logarithms in our code. 

6. Common segment exchange 

There is another possible exchange type, which in- 
volves swapping the allocation of kinks between paths 
for T > Tcs for any common segment at tcs, i e. kinks 
at times t > tcs associated with path A are reassigned 
to path B and those assigned to path B are reassigned 
to path A (common segment means that the two paths 
occupy the same lattice site at the same imaginary time) . 
Common segment exchange is not applicable for all path 
configurations, but has potential to speed up computa- 
tion. In our algorithm, we attempt a common segment 



exchange with probability pcg every time an exchange is 
attempted. Once the attempt is started, we test for the 
existence of a common segment. If there are no shared 
segments, then the update is rejected and no further ex- 
change is attempted on that time step. Otherwise, the 
common segment exchange is accepted with probability 
P(C ^D)= min {1, exp(74(£)) - A{C))]. 



7. Path shift 

To encourage rapid exploration of the path configura- 
tions, it is also convenient to introduce another update 
that is attempted with a small independent probability 
relative to the other updates. In this update, one of the 
paths is shifted by a random number of lattice spacings. 
These updates are useful when the bipolaron is only just 
bound. The update proceeds as follows: 

Attempt the update with an independent probabil- 
ity on each step. Test the configuration of the paths. 
If the paths are in an exchanged configuration, then 
the update is rejected. Otherwise, select one of the 
paths with equal probability. Choose a displacement 
for the path at random, and choose whether to move 
the path left or right with equal probability (impor- 
tant for the detailed balance so the move is reversible). 
The move is accepted with probability P{C —>£)) = 
min{l,exp(A(£)) - A(C))}. 



